In silico Identification of Hypoxic Signature followed by reverse transcription-quantitative PCR Validation in Cancer Cell Lines

Background: Hypoxic tumor microenvironment is one of the important impediments for conventional cancer therapy. This study aimed to computationally identify hypoxia-related mRNA signatures in nine hypoxic-conditioned cancer cell lines and investigate their role during hypoxia. Methods: Nine RNA-Seq expression data sets were retrieved from the Gene Expression Omnibus database. DEGs were identified in each cancer cell line. Then 23 common DEGs were selected by comparing the gene lists across the nine cancer cell lines. qRT-PCR was performed to validate the identified DEGs. Results: By comparing the data sets, GAPDH, LRP1, ALDOA, EFEMP2, PLOD2, CA9, EGLN3, HK, PDK1, KDM3A, UBC, and P4HA1 were identified as hub genes. In addition, miR-335-5p, miR-122-5p, miR-6807-5p, miR-1915-3p, miR-6764-5p, miR-92-3p, miR-23b-3p, miR-615-3p, miR-124-3p, miR-484, and miR-455-3p were determined as common miRNAs. Four DEGs were selected for mRNA expression validation in cancer cells under normoxic and hypoxic conditions with qRT-PCR. The results also showed that the expression levels determined by qRT-PCR were consistent with RNA-Seq data. Conclusion: The identified PPI network of common DEGs could serve as potential hypoxia biomarkers and might be helpful for improving therapeutic strategies.


INTRODUCTION
gene expression signature is a single or combined group of genes whose expression responds to a particular signal or changes in cellular status in a predictable way. Gene signatures are frequently extracted from a set of DEGs by comparing two groups, such as cell lines under different treatment conditions. Gene expression signatures can therefore be used as surrogate markers to comprehend the complexity of pathway activation.
Oxygen deprivation occurs in almost all solid tumors. A shortage of oxygen is the consequence of inadequate oxygen delivery via inefficient tumor vasculature [1] . Hypoxia affects tumor behavior and facilitates tumor progression and metastasis, leading to resistance to conventional chemo-and radiotherapy [2] . Therefore, identifying the key genes regulating cancer cell behavior during hypoxia is essential for developing anticancer agents that efficiently kill tumor cells under hypoxic conditions.
A growing number of studies have identified DEGs A 24 Iran. Biomed. J. 27 (1): [23][24][25][26][27][28][29][30][31][32][33] during hypoxia in different cancer cell types using RNA-Seq analysis [3][4][5] . However, their findings only represent the genetic characteristics of specific tumor cells during hypoxia. In this study, we used RNA-Seq datasets of nine different hypoxic-conditioned cancer cell lines to find hypoxia-related mRNA signatures. Since human cancer cell lines are widely used for better understanding of cancer biology, cancer cell characterization, and anticancer drug discovery [6] , we selected the available RNA-Seq datasets of cell lines to explore the effect of hypoxia on gene expression profiles.
MiRNAs play a central role in regulating gene expression [7] . Kulshreshtha and colleagues [8] described a functional link between hypoxia and miRNA expression. They indicated that miRNAs profile are regulated by hypoxia in a variety of cell types, and their dysregulation is associated with many cancers, making their signature a potential prognostic biomarker [9] . In the present study, common DEGs along with their hub genes among the nine different cancer cell lines were screened during hypoxia. Then we investigated a PPI network and predicted a miRNA-targeted gene network, which might provide a basis for further studies. Our aim was to discover the molecular mechanism underlying the effect of hypoxia and provide potential prognostic markers.

Raw biological data and differential RNA expression analysis
Raw RNA-Seq data of nine hypoxia-conditioned cancer cell lines were retrieved from the Sequence Read Archive (www.ncbi.nlm.nih.gov/geo). Among these datasets, GSE131378 contained four samples of hypoxic-conditioned and four samples of normoxicconditioned A549 cells, while GSE72437 consisted of five samples of hypoxic-conditioned and five samples of normoxic-conditioned BeWo cells. Moreover, GSE78025, GSE81513, GSE84167, GSE13967, GSE149132, and GSE160491 contained three samples of hypoxic-conditioned and also three samples of normoxic-conditioned U78-MG, HCT116, MCF-7, ASPC-1, T47D, and BCPAP, respectively. GSE131379 also comprised of two samples of hypoxic-conditioned and three samples of normoxic-conditioned Hela cells. SAMtools was used to extract raw sequencing reads. The read quality was examined using FastQC version 0.11.2, and low quality bases and adaptor sequences were removed using Trimmomatic version 0.32; the expression level of each transcript was then quantified in transcripts per million using Kallisto [10] . The counts were imported into software R v. 3.4.0 using the tximport R package v. 1.4.0, and the DEGs were identified with a | log2 fold change | ≥1 and a false discovery rate <0.05 using the DESeq2 package in R v. 3.2.3. The UpSetR package in R was employed to find common genes between different datasets [11] . The default values were employed for all the packages.

Construction of a PPI network
Interactions between the common DEGs and other proteins would be useful to fully understand their biological roles. In this study, 23 common DEG PPI network were constructed by Retrieval of Interacting Genes (STRING; https://string-db.org/). Moreover, 23 common DEGs were integrated into the International Molecular Exchange Consortium database (https:// www.imexconsortium.org/) to identify the hub genes information in PPI network [16] . The protein interaction network was visualized using NetworkAnalyst (https://www.networkanalyst.ca) and Cytoscape (3.9.1) [16] . To evaluate the nodes in the PPI network, we adopted several topological measures, including degree (k), MCC, BC, and CC. Since degree (k), BC, and MCC are often used for detecting the hub in a network [17][18][19] , we determined hub genes based on connectivity degree (number of interactions) >10, MCC, and BC using Cytohubba on Cytoscape.

MiRNA interactions analysis
To identify the miRNA-mRNA target interactions, miRTarBase [20] and TarBase [21] (both version 8.0) were employed to collect the miRNA-gene interaction data. Topological analysis based on degree and betweenness centrality as key topological parameters was performed utilizing NetworkAnalyst.

Cell culture for qRT-PCR validation
To validate our findings, we selected four hub genes, including GAPDH, LRP1, ALDOA, and PLOD2 to determine their expression in cancer cell lines (A549, U78-MG, HCT116, Hela, and MCF-7) under hypoxic or normoxic conditions. Cells were purchased from the National Cell Bank of Iran (Pasteur Institute, Tehran, Iran). Cells used in the experiment were cultured in DMEM supplemented with 10% FBS and incubated in a humidified incubator with 5% CO2 at 37 °C. Cancer cell adaptation to hypoxia Cells were seeded in a T25 flask and cultured in DMEM medium supplemented with 10% FBS. The cells were repeatedly incubated in hypoxic conditions in an Anoxomat chamber (Mart Microbiology, Lichtenvoorde, The Netherland; 1% O 2 ) for 4 h and then incubated in a standard culture environment (5% CO 2 and 95% air) at 37 °C for 48-72 h. Cells were treated twice weekly, and hypoxic-conditioned cell lines were generated after 20 exposures to hypoxia [22] .

RNA isolation and qRT-PCR
Trizol reagent (TaKara, Kusatsu, Shiga, Japan) was used for RNA isolation from the cells during normoxia and hypoxia. RNA samples were reversely transcribed to complementary DNA by the QIAGEN Reverse Transcription Kit (Qiagen, Germany). Subsequently, the quantification of cDNA was performed by the qRT-PCR method using SYBR Green Master Mix (Amplicon). The reaction conditions were conducted at 95 °C for 10 min, 40 cycles of 95 °C for 10 s, 60 °C for 30 s, and 72 °C for 30 s. The RPLP0 was used as an internal reference control [23] . Gene expression levels were calculated based on the Delta-Delta Ct relative quantification.

Statistical analysis
Statistical analyses were performed using the student's t-test with GraphPad Prism 8 software (GraphPad Prism, San Diego, CA, USA). The p value was considered statistically significant when it was less than 0.05.

PPI network construction and hub gene selection
Using the STRING database, a PPI network obtained from 23 common DEGs, which was composed of 22 nodes and 25 edges, was constructed and visualized in Cytoscape ( Supplementary Fig. 1). In order to screen the PPI network's interactions with other proteins, which provide important clues about their functions, the PPI network was integrated into the International Molecular Exchange Consortium database. A PPI network composed of 448 nodes and 531 edges was obtained (Fig. 2). Twelve hub proteins, including GAPDH, LRP1, ALDOA, EFEMP2, PLOD2, CA9, EGLN3, HK, PDK1, KDM3A, UBC, and P4HA1, were identified in this network based on degrees (>10), MCC, and BC ( Fig. 3 and Table 3).

Quantitative real-time PCR for DEGs
In order to validate the DEGs identified by RNA-seq analysis, four hub genes, including GAPDH, LRP1, ALDOA, and PLOD2, were selected for analysis via qRT-PCR under normoxic and hypoxic conditions. Primers were designed based on available sequences to amplify the specific altered genes. Primer sequences are shown in Table 4. Based on the qRT-PCR results, the candidate genes were upregulated in A549, U78-MG, HCT116, Hela, and MCF-7 cells under hypoxic conditions (Fig. 5). The expression profiles of four genes confirmed the original transcriptome data obtained by RNA-Seq.

DISCUSSION
Because hypoxic cells are likely to be resistant to chemo-and radiotherapy, it is of high importance to identify the key hypoxia-inducible genes and resistance mechanisms for efficient therapeutic intervention. Moreover, it is well established that miRNA plays a central role in regulating the various biological pathways [24] . Therefore, exploring the role and impact of mRNA and miRNA in cancer cells, especially during hypoxia, could be helpful in cancer diagnosis and treatment.
In the current study, we conducted bioinformatics analysis to identify the candidate key genes and biological pathways among nine different cancer cell lines exposed to hypoxic conditions. Data was extracted from GSE131378, GSE72437, GSE78025, GSE81513, GSE131379, GSE84167, GSE13967, GSE149132, and GSE160491 datasets, among which 23 common DEGs were screened. To our surprise, all the common DEGs were upregulated in all the nine hypoxic-conditioned cancer cell lines. In order to gain some insight into how hypoxia affects the expression of genes at the molecular level, GO and KEGG pathway enrichment analyses were carried out [13,14] . Functional enrichment analysis revealed that the hexose metabolic process, response to hypoxia, and glucose metabolic process were significantly changed. According to KEGG enrichment analysis, 23 common genes were enriched in the HIF-1 signaling pathway, including fructose and mannose metabolism, glycolysis/gluconeogenesis, carbon metabolism, cholesterol metabolism, central carbon metabolism in cancer, and biosynthesis of amino acids. Since it is believed that proteins with more interactions have higher chances of being involved in the essential PPI [25] , the PPI network was constructed and GAPDH, LRP1, ALDOA, EFEMP2, PLOD2, CA9, EGLN3, HK, and PDK1 were identified as the hub genes.
To support our findings, we selected four hub genes (GAPDH, LRP1, ALDOA, and PLOD2) for qRT-PCR validation in A549, U78-MG, HCT116, Hela, and MCF-7 cells under normoxic and hypoxic conditions. Expression patterns of four genes generated by qRT-PCR were consistent with RNA-seq data. Consistently, several studies have found that hypoxia-related genes such as GAPDH, LRP1, ALDOA, EFEMP2, PLOD2, CA9, EGLN3, HK, and PDK1 are upregulated during hypoxia [26][27][28] .    GAPDH and ALDOA are involved in glycolysis. It is widely believed that the overexpression of glycolytic enzymes in a large number of tumors compensates for the increased energy demands and supports rapid tumor growth [29] . However, many glycolytic enzymes have non-glycolytic functions, as well [30] . For instance, overexpressed GAPDH could inhibit caspaseindependent cell death by inducing Bcl-xL upregulation, leading to cancer cell survival and resistance to chemotherapeutic agents [31,32] . Moreover, GAPDH protects cancer cells against chemotherapy by directly binding to the telomeric DNA and prevents the rapid degradation of telomeres [33] . More importantly, GAPDH, which is perceived as a common reference gene, is upregulates under hypoxic conditions. Therefore, using GAPDH as a housekeeping gene should be avoided due to its unstable expression level during hypoxia.
ALDOA and PDK1 are glycolytic enzymes that contribute to the progress of cancer and metastasis [34][35][36][37] . ALDOA overexpression could suppress the expression of proteins responsible for cellcell adhesion and induce the expression of epithelialmesenchymal transition [34] . Chang et al. [34] have demonstrated a feedback loop between ALDOA and HIF-1, by which ALDOA activates HIF-1α/MMP9 and promotes cancer cell invasion. Under hypoxic conditions, PDK1 attenuates mitochondrial respiration and ROS production by inactivating the pyruvate dehydrogenase [38] . Additionally, Gibadulinova et al. [39] have indicated that carbonic anhydrase IX promotes metabolic adaptation to hypoxia through the regulation of PDK1. A number of studies have also revealed that PDK1 overexpression promotes cancer cell metastasis, but the molecular mechanism is unclear [36,37] . Siu et al. [37] have explained that PDK1 expression is associated with ovarian cancer metastasis through the activation of JNK/IL-8 signaling. It has also been displayed that procollagen-lysine, 2-oxoglutarate, PLOD2 promotes migration and invasion of cancer cells during hypoxia. PLOD2, a regulator of collagen cross-linking, is located in the upstream of HK2 and can regulate HK2 expression through the activation of signal transducer and activator of transcription 3 (STAT3) [40] .
In summary, the present study identified hypoxiarelated gene signatures among the hypoxia-conditioned cancer cell lines using RNA-Seq. Our analysis revealed the common hub genes and key pathways in cancer cells under hypoxic conditions. Moreover, we predicted a miRNA signature, among which miR-335-5p had the highest betweenness centrality during hypoxia. To our knowledge, for the first time, our results demonstrate that miR-6807-5p and miR-6764-5p are dysregulated under hypoxic conditions. However, further molecular biological experiments are required to confirm the function of the identified miRNA associated with hypoxia. The results of the present study may provide future directions in identifying the presence of cancer and determining the characteristics of cancer. For instance, hypoxia is a characteristic feature of cancer, and the hypoxia signature identified in this study, as well as predicted miRNAs might be helpful to detect the hypoxic state of cancer cells. Hypoxia is common in majority of malignant tumors and an attractive therapeutic target. As hypoxia targeted treatment are effective in patients with the most hypoxic tumors, hypoxic signature might be useful for developing proper treatment, such as engineered oncolytic viruses that could be utilized to control or regulate the biological interactions responsible for the functioning or malfunctioning of cancer cells during hypoxia [22] .

Ethical statement
Not applicable.

Data availability
The raw data supporting the conclusions of this article are available from the authors upon request.

Author contributions
SS: conceived, designed the analysis and performed the analysis; GB: performed bioinformatics analyses; AA: drafted or provided critical revision of the article; KA: conceived and designed the study, supervised the data analysis and interpretation. All authors have read and approved the final manuscript.